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A recent technique, proposed to alleviate the "sign problem disease", is 
discussed in details. As well known the ground state of a given Hamiltonian H 
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0> , can be obtained by applying the imaginary time propagator e to a given 



trial state ipT for large imaginary time r and sampling statistically the propa- 
gated state TpT = e~^'^'ipT- However the so called "sign problem" may appear 



O . in the simulation and such statistical propagation would be practically impos- 

sible without employing some approximation such as the well known "fixed 
node" approximation (FN). This method allows to improve the FN dynamic 
with a systematic correction scheme. This is possible by the simple require- 
ment that, after a short imaginary time propagation via the FN dynamic, 
a number p of correlation functions can be further constrained to be exact 
by small perturbation of the FN propagated state, which is free of the sign 
problem. By iterating this scheme the Monte Carlo average sign, which is 
almost zero when there is sign problem, remains stable and finite even for 
large r. The proposed algorithm is tested against the exact diagonalization 
results available on finite lattice. It is also shown in few test cases that the 
dependence of the results upon the few parameters entering the stochastic 



technique can be very easily controlled, unless for exceptional cases. 
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I. INTRODUCTION 

In the last few years an enormous progress in the computational techniques has also been 
accompanied by better and better performances of modern computers. All these develop- 
ments have certainly contributed to determine a "feeling" that the many body problem of 
solving a strongly correlated Hamiltonian, with many electrons on a reasonably large system 
size, is becoming possible with some computational effort. 

The various numerical methods, like e.g. to find the ground state of a physically interest- 
ing Hamiltonian, can be classified in two main branches developing from two root methods: 
the exact diagonalization technique (ED) and the variational Monte Carlo method (VMC). 

The first technique is a brute force diagonalization of the Hamiltonian matrix, which 
represents a prohibitive task for large number of electrons as the linear dimension of this 
matrix grows exponentially with the number of electrons and the system size. The use of 
spatial symmetries and the very efficient Lanczos technique have made recently possible 
the exact ground state evaluation of up to ~ 30 electrons on simple lattice Hamiltonians 
like: the Heisenberg model [|l|, the t — J model 0, the Hubbard model and similar ones 
01 . However this system size is far from being enough for the determination of the physical 
thermodynamic properties of the various models. 

A remarkable development of the ED like methods, is certainly the so called density 
matrix renormalization group technique (DMRG). In this case the ground state of a huge 
Hilbert space Hamiltonian is sampled by a small basis set that is iteratively improved by 
using the renormalization group idea. In one dimension this technique allows to have for 
instance the numerically exact solution of the Heisenberg spin S = 1 model for the infinite 
size m. Recently DMRG has also been extended for high accuracy calculations on simple 
molecules 0. 

The second root of development starts from the VMC technique 0. The VMC allows 
to sample statistically a variational wavefunction ipcix)-, defined on a given basis set, whose 
elements {%} are represented by simple configurations defined typically by the electron 



positions and spins. In the most simple formulation the VMC sampling can be obtained 
by accepting a new configuration Xn+i from a given one Xn if a random number ^ between 
zero and one satisfies C, < lipcixn+i) /ipGixn)]"^ , otherwise Xn+i = Xn- This simple Metropolis 
algorithm generates states x„ that, after some equilibration, are distributed statistically 
according to the square of the variational wavefunction. Then physical expectation values 
of operators O'' - such as pair correlations, electron number, total spin square, energy etc. 

can be easily obtained on the given variational wavefunction, provided the local estimator 
O^ = /^ u) of the correlation function O^ can be computed in an efficient way. This is 
typically the case since the configuration basis is particularly simple so that {iPg\x) can be 
easily computed, and also 0''\x) can be expanded in a few (less than the square electron 
number) configurations, for one and two body correlations. 

The iterative rule determining a new configuration Xn+i starting from a previous one x„, 
and depending also on a random number, defines a Markov chain which allows to obtain 
statistical estimates of the above expectation values. This is possible even if the dimension 
of the Hilbert space is very large, such property representing the most important advantage 
of the statistical methods. 

From this point of view the Green function Monte Carlo (GFMC) technique can be 
considered a development of the VMC because it allows to sample statistically the exact 
ground state of a many body Hamiltonian H, instead of being restricted to the variational 
wavefunction. In the GFMC the ground state is statistically sampled by a set of M walkers 
(wj,Xi), i = 1,---,M, i.e., at each configuration Xi is associated a weight Wi in order to 
represent a simple element Wi Xi of the large (or even infinite) Hilbert space. In this case 
a Markov chain, which is slightly more complicated than the variational one, can be easily 
defined. As it will be shown later the new configurations and weights (wj,Xj)„+i depend 
only on the previous weights and configurations {wi,Xi)n and M random numbers C,i- This 



iteration is equivalent statistically to a matrix-vector product 

^n+l{x') = Y^ G^,^^^n{x) (2) 

X 

where Gx\x is the lattice Green function simply related to the Hamiltonian matrix elements 
in the given basis 

^ x' X ~~ ^^^x' X — ^x' X \y) 

and A is a suitable constant, allowing the convergence of @) to the ground state of H for 
large n. At each Markov iteration n the state ipn{x) is sampled statistically by the walkers, 
which may be even a large number, but typically a neglectable fraction of the total Hilbert 
space dimension. 

In the statistical iteration the weights Wi of the walkers increase or decrease exponentially 
so that after a few iterations most of the walkers have an irrelevant weight w and some kind of 
reconfiguration becomes necessary to avoid large statistical errors. The process to eliminate 
the irrelevant walkers from the statistical sampling is called "branching" . This amounts for 
instance to duplicate a walker with large Wi in two walkers with half the weights Wi/2 acting 
on the same configuration, or drop the walkers with too small weights. If properly done 
this kind of process does not introduce any bias but the number of walkers is not constant 
during the corresponding Markov chain. For practical purposes it is necessary therefore 
to control the walker population number otherwise the simulation exceeds the maximum 
available memory or terminates for lack of walkers. This statistical reconfiguration instead 
introduces some amount of bias. Recently a rigorous and simple way to work at finite 
number of walkers has been proposed, which simplifies the GFMC technique by controlling 
and eventually eliminating the bias due to the finite number of walkers |p. 

With a slight generalization of the previous simple technique it is also possible |10| to 



alleviate the "unfamous sign problem" , which occurs when the matrix elements of the lattice 
Green function G^^x are not always positive definite. In this case the iteration (|^) can still 
have a statistical meaning at the price that the weights Wi of the walkers are no longer 



restricted to be positive. It then happens that the average sign (s)„ = -'^ at a Markov 
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iteration n is exponentially decreasing with n, implying a dramatic decrease of the signal 
to noise ratio for all correlation functions. A remarkable improvement of the GFMC on 
a lattice was the extension of the fixed nodes (FN) approximation to lattice Hamiltonians 
0. In this case the "dangerous" negative off-diagonal elements of the Green function are 
neglected and stable simulations with always positive walker weights Wi can be performed 
at the price of obtaining an approximate solution of the ground state wavefunction. 

The Green function Monte Carlo with Stochastic Reconfiguration (GFMCSR) ||10[ rep- 
resents a successful attempt to improve the FN, with a stable simulation without any sign 
problem instability. In this scheme, better and better approximations of the ground state 
correlation functions may be obtained, by performing controlled Markov chain simulations 
with average walker sign (s)„ very close to 1 for each iteration n. For the sake of simplicity 
we restrict the forthcoming derivation to lattice Hamiltonians but the basic ideas can be 
straightforwardly extended to the continuous case. This method is based upon the simple 
requirement that after a few iterations of (0) via the approximate FN dynamic, a number 
p of correlation functions can be further constrained to be exact by properly small pertur- 
bations of the propagated FN state ^^■^■^, which is free of the sign problem. By iterating 
this process the average sign remains stable even for large n and , in this limit, the method 
has the important property to be in principle exact if all possible correlation functions are 
included in this correction scheme of the FN. 

In the first five sections we review the basic steps of the GFMC for the general case 
when the sign problem affects the practical implementation of the algorithm. In Sec. |V| we 
introduce the Stochastic Reconfiguration (SR) idea and in Sec. |Vll|we prove the fundamental 



theorem, which justify the approximations used to get rid of the sign problem. In the 
remaining sections we present the details of the algorithm and some test results, useful 
to understand how to implement the numerical algorithm, for an efficient and controlled 



improvement of the FN, even for large system sizes. 

II. THE GFMC TECHNIQUE 

From a general point of view the ground state tpQ of a lattice Hamiltonian H can be 
obtained by iterating the well known power method (0) so that ^„ —* tpQ for large n, 
provided the initial state ipT at the first iteration of Eq. (^ {ipn = ipT for n = 1) is a trial 
state non-orthogonal to the ground state ipQ. 

A stochastic approach is possible if one can sample statistically the matrix-vector iter- 
ations d^). This is particularly important since for large systems only few power iterations 
can be applied exactly in the most fortunate cases. The important property that allows 
a statistical approach is that physical lattice Hamiltonians are represented by very sparse 
matrices. Though the total number of non-zero elements of G^'.^ is prohibitive, the number 
of non- vanishing entries in each column is a neglectable fraction - of the order of the electron 
number - of the total Hilbert space dimension. Thus all the non-zero G^^x for fixed column 
index x can be computed even for large size. 

It is therefore natural to define a basic element of the stochastic approach: the so called 
walker. A walker is determined by an index x corresponding to a given element \x) of the 
chosen basis and a weight w. With a stochastic approach the walker "walks" in the Hilbert 
space of the matrix H and assumes a configuration wx according to a given probability 
distribution P{w,x). 

The task of the GFMC approach is to define a Markov chain, yielding a probability 
distribution Pn{w,x) for the walker which determines the iterated wavef unction ipn'- 

V'„(x) = {x\^n) = / dwwPn{w,x) . (4) 



III. IMPORTANCE SAMPLING 

One of the most important advantages of the GFMC technique is the possibihty to reduce 
the variance of the energy by exploiting some information of the ground state wavefunction, 
known a priori on physical grounds. In order to understand how, we simply note that the 



power method is not restricted to symmetric matrices. Following Ceperley and Kalos ||Tl 
one can consider in the iteration (@) not the original matrix G, but the slightly more involved 
non-symmetric one 

Gx',x = '4'g{x')Gx',x/'^Pg{x) , (5) 

where ipc is the so called guiding wavefunction, that has to be as simple as possible to 
be efficiently implemented in the calculation of the matrix elements and, as we will see, as 
close as possible to the ground state of H. Here and in the following we assume that the 
guiding wavefunction is always non- vanishing for all x. It is obvious that G, though being a 
non-symmetric matrix, has the same spectrum of G as to any eigenvector ipk{x) of G with 
energy A — E^ corresponds a right eigenvector of G equal to 'ipG{x)'^k{x) with the same 
eigenvalue. 

As shown later on, by sampling statistically the iteration (^ with G instead of G the 
walkers {w,x) will be distributed for large n according to iPo{x)4'g{x), namely '?/'„(x) ex 
ipo{x)ipG{x) in Eq. (||). In order to evaluate the ground state energy, it is then enough to 
average the so called local energy. 

Ex = iM^ = Y: Mx')Hx'JMx) , (6) 

over the statistically sampled walkers, because obviously: 

J2x Exll>oix)^Gix) 



{E.) 



^"'^^ ExM^)M^) 
{Mi^G) 

Thus if ifjc is exactly equal to the ground state of H, by definition Ex = Eq, independent 
of X, as {iPg\H = Eq{iPg\ in ®- This is the so called zero variance property satisfied 
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by the method. Namely if the guiding wavefunction approaches an exact eigenstate of H, 
the method is free of statistical fluctuations. Of course such a fortunate situation is not 
common at all, but by simply improving the guiding wavefunction the statistical fluctua- 
tions of the energy are much reduced, leading to more efficient simulations. This property, 
rather obvious, is very important and non-trivial. Many methods, in fact, such as the path 
integral Monte Carlo, suffer of statistical fluctuations even if an exact information of the 
desired eigenstate is known. For Hamiltonians affected by the sign problem it is particularly 
important to work with a method which depends strongly on the quality of the initial guess 
of the ground state represented in the GFMC by the guiding wavefunction. This helps a 
lot because by the simple and successful "trial and error" strategy one can systematically 
improve the guiding wavefunction and gain information about the ground state. 

In general after the transformation (|]) all mixed average correlation functions defined 
by linear operators O^ 

(^gI^o) ^ ^ 

are easily accessible by GFMC. The local estimator corresponding to Eq. (|^) is, analogously 
to (D, given by 

0.' = EV'g(x')OV^g(x), (8) 

x' 

exactly as in the variational approach ([l|). This expression represents just the sum over all 
the possible matrix elements connected to x of the transformed operators O^ with matrix 
elements 

0%,. = i^G{x')Ol.J^G{x) , (9) 

namely O^ = J2x' O^' x- I^ order to implement the "importance sampling" strategy it is 
sufficient therefore to replace all the matrices involved O^, ^ including the Green function 
G with the transformed ones O'' and G (^, and in all previous expressions the guiding 
wavefunction disappears. Thus the method can be considered a general method to find the 
maximum eigenvalue and eigenvector of a generic (non-symmetric) matrix G. 

9 



In the following, for simplicity of notations, we put a bar over the symbols corresponding 
to all the transformed matrices (|]) and (H). 

IV. SINGLE WALKER FORMULATION 

In general the distribution P„(w,x) is sampled by a finite number M of walkers. Let us 
first consider the simpler case M = 1. In order to define a statistical implementation of the 
matrix multiplication (|^), the standard approach is first to determine the Green function 
matrix elements Gxi,x connected to x which are different from zero. These matrix elements 
can be generally written in terms of three factors 

Gx',x = Sx',xPx\xbx (10) 

where bx is a positive normalization factor, Sx'^x takes into account the signs of the Green 
function and Px'^x is a stochastic matrix. All these terms will be defined explicitly below. 

The basic step of the GFMC method on a lattice is to define properly the matrix Px',x, be- 
cause it represents the only term in the decomposition (|TDp that allows to select statistically 
only one configuration among all the possible ones x' connected to the single configuration 
X of the walker by the Green function application (^. Therefore Px\x has to represent a 
probability and is restricted to be i) normalized J2x' Px'.x = 1 and ii) with all positive matrix 
elements Px',x > 0. This is just the definition of a stochastic matrix (see Appendix ^. Since 
the matrix elements of G are not restricted to be positive (sign problem) px'^x is more clearly 
defined in terms of an appropriate Green function G^-^-^ with all positive matrix elements. 
Even if the latter restriction may appear rather strong, it is however possible that for large 
n the approximate propagation of the state ip^^-^ by the Green function G'^-^-'^ is not far from 
the true propagation of ipn by the exact Green function G in Eq. (0). G^.{{, needs not to be 
normalized, as its normalization can be included in the definition of the positive constant 

bx = T.GZ (11) 

x' 

so that 
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gZ = P^',^K . (12) 

The typical choice for G*^-^-^ is given by the absolute value of the matrix elements of G, 
Gx'x ~ l^^'.^l' ^^^ ^^^^ ^^ ^°^ ^^^ optimal choice as it will be discussed below. 

Since the most stable right eigenvector ipQ^(x) of a positive definite Green function - 
like G^^^ - can be chosen positive -ipQ (x) > 0, it is important to implement importance 
sampling by the transformation (^) with a guiding wavefunction with signs as similar as 
possible to the ones of the ground state of H, so that the Green function G has its most 
stable right eigenvector iljG{x)'ipo{x) > for most configurations x. In this case there are 
good chances that the latter state is well approximated by the positive vector tp^^^ [x) > 0, 
generated by G'^-^^ for large n. In order to fulfill better the latter requirement, here we follow 
a recent development of the FN on a lattice, and we choose for G'^^^^ the FN Green function 
(with importance sampling): 

GZ = A5..,. - ma ■ (13) 

The constant shift A has to be large enough that all the diagonal elements of G'^-^'^ are strictly 
positive. This is possible in general for the diagonal elements. The full Green function G^^-^ 
is defined in a way that the ground state of the Hamiltonian H'^^^, is a variational state of 
H with an energy better than the guiding wavefunction one [|l^. Contrary to the standard 
FN method, that neglects all the matrix elements of H that cross the nodes of the guiding 
wavefunction, namely the ones with H^i^x > 0, we adopt here a slight modification of H^^^ 
defined with non-zero matrix elements (but with opposite sign) when H has the positive 
ones. The generalization of the above "FN theorem" to this case is straightforward and is 
reported in the Appendix ^. 

The appropriate matrix elements of H'^^^ are obtained by reversing the sign of the positive 
off-diagonal matrix elements of H and by multiplying them by a constant 7 > 



/It-' X if -fir' r "^ U 



'x',a; 



H:ji = { ' - (14) 

—'jHx^x if Hx\x > 
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and the diagonal ones are 

i7:f/ = i7.,. + (l + 7)Vsf(a;), (15) 

where the diagonal sign-Hip contribution is given by [ p^ : 

Vsf(x) = Y. H^',- ■ (16) 

H^i ^>0 and x'j^x 

Notice that there is no difference between the diagonal elements of the Hamiltonian H'^^'f 
(H) and the ones of the transformed matrix ff*^-^-^ (H), as defined by Eq. (^. 
The equality (|1^) holds if the factor s^'^x is given by: 



^x': 



1 


if G^'.x > 


-1/7 


if G;.,^, < . 


A-Hx.x 
A-H^x^J 


if x' = X 



(17) 



The appropriate stochastic process relative to the Hamiltonian H can be defined in the 
following three steps, simply by allowing the weight w of the walker to be also negative: 

1. Given the walker {w,x), change the weight by scaling it with bx'- 

w -^ hxW . 

2. Generate randomly a new configuration x' according to the stochastic matrix Px\x- 

3. Finally multiply the weight of the walker by s^'^x'- 

Without the latter step, one is actually sampling the Hamiltonian H'^^^ , which we expect 
(or assume) to have a ground state close to the one of H, for suitably chosen guiding 
wavefunction. During the Markov iteration (MI) it is straightforward therefore to update 
both the weight w associated to the true Hamiltonian and the one w^-^-'^ associated to the 
approximate one W^^f . From now on the walker will be therefore characterized by the triad: 

{w,w'^^-^ , x) . 
12 



The previous MI allows to define the evolution of the probability density to have the 
walker with weights w and w'^^^ > in the configuration x, namely: 

p„,,K ^.^//^ x') = Ei|r^^"(r^,^,-) • (18) 

X ^x\^x' ,x\ ^x^x' ,x "x 

The first momentum of the distribution P over w gives information about the state ipni^) 
propagated with the exact Green function G and the state ip'^^^{x) propagated with the FN 
Green function G'^^^ , namely: 

tpn{x) = I dw^ff f dwwPn{w,W^ff,x) , (19) 

^^//(x) = f dw^ff ( dww^ff Pniw.w^^f.x) . (20) 



In fact it can be readily verified using ([T8|) that the above expressions for ■?/'„ and '4)^^^ i 
satisfy the iteration condition (^ with G and G^^^ respectively. 

At this stage the algorithm is exact and the MI allows to sample the ground state 
of H (with sign problem) and H'^^^ (with no sign problem) within statistical errors, that 
unfortunately may be very large, and increasing with the iteration number n, especially 
when there is sign problem. 

In order to have an idea on the origin of the sign problem let us discuss the follow- 
ing example. Suppose that H'^! {. = —\H\x',x for the off-diagonal elements and H'^-^^ = H 
otherwise. The only information of the difference between the matrix H with respect to 
jj^ff is given by the sampling of the sign. In particular it is easy to realize that in 
this case w'^-^-^ = \w\ at each Markov iteration n. Then at a given iteration n we get 
/ dw^^f f dww^^f Pniuj, w'^f^, x) = f dw'^f^ f dw \w\ Pn{w, w^-f^ , x) ~ (A — Eq-^)", where -Eg 
is the ground state energy of JP^^-^ which is obviously below the ground state energy Eq of 
H. We obtain therefore the basic instability related to this Markov process, known as the 
sign problem, which, as well known, is particularly difficult for fermion systems: 



ExJdw^^f JdwwPn{w,w^ff,x) , A-^i 







(^n) = ^v.......;.......,;v' ....;:^ - i^^^r- (21) 



J2xl dWff Jdw\w\Pn{'w,w''ff,x) A-E, 







The latter relation shows that, for large n, walkers with positive weight w > cancel almost 
exactly the contribution of the walkers with negative weight w < leaving an exponentially 
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smaller quantity which is obviously difficult to sample. In this case only few power iterations 
ra ~ 10 are possible |[T^ and for large system size this is by far not sufficient even for a minor 
improvement of the initial guess ipc- It is important to emphasize that this instability does 
not even depend on the guiding wavefunction because the latter cannot change the spectrum 
of H and H^^^ defined above. 

By iterating several times the MI even for a single walker, the resulting configuration 
(w7, x) will be distributed according to the ground state of H and by sampling a large number 
of independent configurations we can evaluate for instance the ground state energy 

where the brackets (. . .) indicate the usual stochastic average, namely averaging over the 
independent configurations. 

The configurations Xn that are generated in the Markov process are distributed after 
many iterations according to the maximum right eigenstate of the matrix Px',x (as, if we 
neglect the weights of the walkers, only the matrix p is effective in the matrix product (^). 
This state is in general different from the state ipG{x)'4'f){^) we are interested in. So after 
many iterations the sampled configurations x„ are distributed according to an approximate 
state and we can consider this state as a trial state ipx for the initial iteration n = 1 in the 
power method (0). At any MI n we can compute the weight of the walker assuming that 
L iterations before its value was simply w = 1. In this way it is simple to compute the 
resulting weight of the walker with L power Green function G applications: 

L 

Ti _L J_ ^n — j ^7i—j-\-l i^n—j ' \ J 

i=i 
Therefore for instance, in order to compute the energy with a single Markov chain of many 
iterations, the following quantity is usually sampled 

^0 - V rL ' ^^^^ 



with L fixed |14| 
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This would conclude the GFMC scheme, if averages over the weight variable G^ were 
possible in a stable and controlled manner. However there are two important drawbacks 
for the single walker formulation. The first one arises because the weight G^ of the walker 
grows exponentially with L - simply as a result of the L independent products in Eq. (^) - 
and can assume very large values, implying diverging variances in the above averages. This 
problem has a very well established solution by generalizing the GFMC to many walkers 
and introducing a scheme that enables to carry out walkers with reasonable values of the 
weights, by dropping the irrelevant walkers with small weights and splitting the ones with 
large weights. Recently a simple formulation of this scheme was defined at fixed number 
of walkers p| in a way that allows to control efficiently the residual bias related to the 
finite walker population, as discussed in the introduction. The second drawback is the more 
difficult one and is due to the unfamous sign problem. The average sign (s^) = yc"' 2| 
vanishes exponentially with L as in Eq. (0). In the formulation of Ref. 0] this problem 
looks quite similar to the first simple one. As we will see later on, some kind of remedy 
can be defined by a simple generalization of the SR which is useful in the case with no sign 
problem. 

V. CARRYING MANY CONFIGURATIONS SIMULTANEOUSLY 

Given M walkers we indicate the corresponding configurations and weights with a couple 
of vectors (w, x), with each vector component [wi^wl ,Xi) i = 1,---,M , corresponding 
to the i^^ walker. Following Q] it is then easy to generalize Eq. (|18|) to many walkers by 
the corresponding probability Pn{ui,x) of having the M walkers with weights and configu- 
rations {w, x) at the iteration n. Similarly to the single walker formulation the propagated 
wavefunctions ipn{x) and ip^^^{x) with the true Green function G and the approximate one 
G'ff read 



(25) 

rj^ix) = J[dw]E^^Hr^^Pn{w,x) 
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where the symbol J[dw\ indicates the 2M multidimensional integral over the {wi^wl ) vari- 
ables i = 1, ■ ■ ■ , M ranging from — oo to oo. Equations (|25D are very important because they 
show that the propagated quantum mechanical states ipn and ip^^^ ) which are sampled sta- 
tistically, do not uniquely determine the walker probability function Pn{w.,x)- In particular, 
it is perfectly possible to define a statistical process, the SR, which changes the probability 
distribution P„ without changing the exact information content, i.e., the mentioned propa- 
gated states ipn and ip!^^^ ■ In this way a linear transformation of P„, described by a simple 
kernel function X{v/,x!;w,x), will be explicitly given: 

Plrim!, x!) = [[dw] J2X{n/, x!; w, x)Pn{w, x) . (26) 

When there is no sign problem it is possible to define the function X [§ in a simple way 
by requiring that the weights w'^ = Wj are all equal to J2j Wj/M after the SR. In this case 
the algorithm is exact, and allows to perform stable simulations by applying the SR each 
few kp iterations. Further, by increasing the number of walkers M, the exponential growth 
in the variance of the weights Wj can be always reduced and systematically controlled. In 
fact for large enough M it is possible to work with L sufficiently large {L oc M) and obtain 
results already converged in the power method iteration (|^) and with small error bars. 

VI. STOCHASTIC RECONFIGURATION, STABILIZATION OF THE SIGN 

PROBLEM 

In order to avoid the sign problem instability, at least in an approximate way, we can 
follow the previous scheme as before by using the following function X that defines the SR 

X(«;',x';m,x) = n ( ^^^f<A 5{w[ - /?-^sgnp.,) Siw^ ' - W^) , (27) 

where /3 = -r^-^^ — ^ is the average sign after the reconfiguration which is supposed to be much 
higher to stabilize the process. The kernel ( pTj) has a particularly simple form since the 
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outcome variables x' and w'a are completely independent for different j values. In particular 
it is possible to integrate easily each of the M factors of the kernel in the variables w', 
Wj and to sum over the configuration a;', the result being one, as it is required by the 
normalization condition for P' in p^). 

After the SR the exact information sampled is obtained by using Eq. (^) with P' instead 
of P. We define the corresponding quantum states '^/'^(a:) and ipn^^'i^)^ ^^^ ^-f^ being exact 
whenever 

i:'^{x)=M^)- (28) 

After the SR the new configurations x'^ are taken randomly among the old ones {xj}, ac- 
cording to the probability ^e^^, defined below in terms of the given weights {wj}, {Wj} 
and configurations {xj}. After that the weights w'^ are changed consistently to (^) in 
w'i = f3~^ Ij ^ sgnpx'. and the FN weights, since are restricted to be positive, are defined 
by taking the absolute value of the previous ones w^ = \w'^\. The coefficient (3-y^j^ 
guarantees the normalization of the two quantum states after and before the reconfigura- 
tion, namely J^xi^'ni^) = J^xi^nix). This coefficient /3 represents also the expected average 
walker sign (s)' = ^ \ ^ after the reconfiguration. It is supposed to be much higher than 
the average sign before the reconfiguration (s) = rr ^\ \ , so that a stable simulation with 

to to to ^ ' l^j \^3\ 

approximately constant average sign (s)' can be obtained by iteratively applying the SR 
each few kp steps of the power method iteration @. 



In the actual implementation of this algorithm (see Sec. |VIII| for the details) the weights 



are reset to unit values after the SR: w[ = sgap^' and w^ = 1, whereas only the overall 
constant /3~^j — , common to all the different walkers, is stored in a sequential file. Then, 
as in the single walker formulation, at any given iteration n, we can assume that L iterations 
before the trial state ipT is given by the equilibrium distribution of walkers with unit weights 
Wj = sgnpxj ■ Therefore in order to obtain the weights predicted by the Eq. ( P^ ) for L power 
method iterations starting from ipT it is enough to multiply the previous L/kp saved factors 
fn = P~^ If ^ ■ This yields a natural extension of the factors G^ (|2^) to the many walker 
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case 

G'n=n/n-fcxfc, (29) 

fc=l 

and the corresponding mixed average correlation functions are obtained by averaging the 
local estimators over all the iterations n just before the SR (i.e. n is a multiple of kp) 

(O-) = y^^i'f-' , (30) 

where, in the above equation, the weights Wj and the local estimators O^ . are evaluated only 
before the SR. 

The only left quantity to define properly the whole algorithm consistently with Eq. (|27D 
are the important coefficients p^ which have not to be assumed positive. These coefficients 
may depend on all the weights Wj, the configurations Xj and the FN weights Wj . 

The choice px = Wj is exact in the sense that ipni^) = V'n(a^), and coincides with the 
one for the case with no sign problem ^. However this choice is obviously not convenient, 
because this reconfiguration will not improve the sign, which will decay exponentially in the 
same way. 

Instead, in the case with sign problem, we can parameterize the coefficients p^ by as- 
suming they are close enough to the positive definite weights {Wj }, the ones obtained with 
the FN Green function G^^^ . The rational of this choice is that, though the weights Wj 
may be occasionally very different from the exact weights Wj - namely the sign can be wrong 
- they sample a state ipn^'^i^) "which is supposed to be quite close to the exact propagated 
state ipn{x). This condition is clearly verified for an appropriate choice of the guiding wave- 
function ipG, which makes the FN accurate. Then we assume that small perturbations over 
the state ip^^^lx) may lead to fulfill the equality ( pHD with an arbitrary small error. In the 
case with sign problem in fact, we release the exact SR condition (p8D to be satisfied within 
some error. This error will affect the equilibrium walker distribution Pn for large n, but 
there will be no problem if this error i) is small and ii) can be reduced within the desired 
accuracy. 



In the most simple and practical formulation we require that the average energy before 
and after the SR coincide 

J2 H^>^^1pnix) = J2 Hx',xi^'n{x) (31) 



(the denominators in the mixed averages (^ are already equal by definition as Z^xV'nl 



X] 



J^xi^nix) for the chosen (3 in (p7|)). Then we define 
and 



E. 






. w ■ 

^3 3 






■'3 



where E^. is the local energy (^) associated to the configuration Xj. Thus E represents the 
estimate of the average energy correctly sampled with the sign, whereas E^ff is the one with 
no sign problem. In order to satisfy the requirement (|3ll) we just determine a by 



A simple calculation shows that with this reconfiguration, that clearly improves the sign, 



where E^jj = ^ ^^^ ,.jy ^ is the average square energy over the positive distribution w 



the value of the energy (the mixed average energy) remains statistically the same before 
and after the SR (see next Section and Appendix [C]). It is clear however that this is not 
enough to guarantee convergence to the exact ground state, because fulfillment of (^TJ) does 
not imply the exact equality (^). We can improve the definition of the constants Px by 
including an arbitrary number p of parameters with p « M 

p,^ = wf{l + a,{Ol^ - Olff) + ■■■ + a,{Ol^ - O^^^)) (34) 

proportional to the fluctuations O^ — O^^r of p different operators O'' with corresponding 
local estimators O^. = ^ ?} , '^■'^ for /c = 1, ■ ■ ■ ,p, and average value over the positive weights: 
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^eff — — y^ ' e// ^ • With the more general form ( p4D for the coefficients px^ it is possible to 
fulfill that all the mixed averages for the chosen p operators - not only the energy - have 
the same value before and after the SR: 

5:0,%^.(x) = E0,^,<(a;). (35) 

In general the reference weights Wj in Eq. (|3^ may be also different from the ones generated 
by the FN Green function, the only restriction is that Wj > for each walker j (see 
Appendix 0). 

It is proven in the next Section that in order to fulfill exactly the SR conditions (|35D it 
is sufficient that the coefficients px are chosen in a way that 



T.P.pl^ i:wpi. 






(36) 



which can be fulfilled with a solution of a simple linear system for the unknown variables 
Ok, for A; = 1, ■ ■ ■ ,p, as described in the Appendix y. The conditions (^) are much simpler 
because they can be satisfied at a given iteration of the Markov process. The theorem, 
proven in the next section, guarantees that the exact (|35|) are implied by the constraints 
(^) after the complete statistical average over the probability walker distribution P„. 

Thus, asymptotically, by adding more and more parameters {aj}, we can achieve 
ip'nix) = ipn{x) strictly, since the distribution ipni^) is completely determined by its cor- 
relation functions. The proof of this important statement is very simple. Consider first 
the diagonal operators. All these operators may be written as linear combinations of the 
"elementary" ones O^?^ = Sx'^xSx,xo acting on a single configuration xq, plus at most some 
constants. If conditions (|35D are satisfied for all the elementary operators it immediately 
follows that ip'^lxo) = ipni^o) for all Xq, which is the exact SR condition (PS]). 

Then it is simple to show that the coefficients Px , determining P^ and ipl^, are invariant 
for any constant shift of the operators O'^. Further with a little algebra it turns out that 
these coefficients Px do not change for any arbitrary linear transformation of the chosen 
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operator set: O^' = J2kLk',kO^ (with real L and det L 7^ 0) (see Appendix P 1|) . Thus the 
proven convergence of the GFMCSR is obtained for any sequence of diagonal operators, that, 
with increasing p, becomes complete. For non-diagonal operators Ox',x note simply that they 
assume the same mixed average values of the equivalent diagonal ones O J°f = 6x',x Y.x' Ox',x- 
Thus the proof that GFMCSR converges in principle to the exact solution is valid in general 
even when non-diagonal operators, such as the Hamiltonian itself for the energy, are included 
in the conditions (|35|) □. 

VII. FORMAL PROOF OF THE GFMCSR CONDITIONS 



As stated before the SR conditions ( p5D read 

EO,%<(x) = EO,^',.^n(x), (37) 



X' ,x 



ioY k = 1,- ■ ■ ,p, with the normalization one J2x V'n(^) = ^x i'nix). 

The wavefunction ip'ni^) after the SR conditions defined by (|27|) can be explicitly written 
in terms of the original walker probability distribution. To this purpose we single out in the 
definition oi 'ip'^{x) 

ij'^ix) = /[rfml E^n(m',x')5^^^, (38) 

x' 

a term k in the above summation over j which gives an additive contribution to ip'^, namely 
< = ]g Efc WJk with 

{i^M}k = [[dm:] E f[dm]J2X{w:,x!-w,x)Pniw,x)Sx,x'wi , (39) 

•^ x' -^ X 

where in the above equation we have substituted the definition of P' in terms of P given by 
Eqs. (|26| ) and (^). In the latter equation it is easy to integrate over all variables WpWj ', x'j 
for j y^ k using that the kernel X is particularly simple as previously discussed. Then, the 
remaining three integrals and summations over w'/^, wl ', x'^ can be easily performed using 



the simple S functions which appear in the kernel X and the definition of /3 = ^ ^, "'•', , so 



J2 P: 

wiiicii appear m me Kernel yv aiiu iiie ueiuiiLion oi /j = ' 
that one easily obtains 
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WA^)}, = [[dw]Y.Pnim,x)^^sgnp^?^2p;^ . (40) 

J X ^'^ Z^j Pxj 

It is important to remark that, in the above equation, the sign of px (sgnp^;) depends only 
on the configuration x chosen among the old configurations Xj, determining the vector x in 
Pn{w,x). In particular if there are more walkers acting on the same configuration {xj = x 
for more than one j) sgnp^ is the same for all the corresponding indices, as implied by the 



eff 



definition (p^ ) of p^^^ and the condition Wj > valid for all j. We can therefore replace in 



y . \px ISx X ■ / - Px -^x X ■ 

general sgnp^; — W^^ — '-^ = — ^-^ — '^- and obtain the closed expression for ipn{x) after the 



simple summation on the index k: 

i^M = lmT.Pn{ni.x){^)^^ . (41) 

J X ^"^ 2-1 j Pxj 

Then the normalization condition J2x^nix)= J[dw] J2x Pn{iil,x){ ^j ^ ) = J^xi^nix) eas- 
ily follows. On the other hand the left hand side of Eqs. (^) can be also computed easily, 
yielding 

j:0'x',x^ni^) = [[dw]j:Pn{w,x){?^)?^^^ , (42) 

x',x ^ ^^^ 2^jP^j 

where O^. = J2x' Ox'^xj is the mixed estimator of the operator O'^. 

Finally, by substituting the conditions ( p6D into the previous equation, one obtains 

y: o.^,.v^;(x) = j [dw] y: Pniw,x) ' ^ ^- = y: ot,xMx) , (43) 

x\x X x' ,x 

which proves the statement at the beginning of this section. 

A. Optimization of the weights 

The definition of the weights px that satisfies the SR conditions (^) is highly arbitrary 
because as we have mentioned before the probabilities P„ and P^ do not uniquely determine 
the quantum states ipn and ip'n that are subject to the conditions (|35|) . In this sense there 
may be different definitions of the weights Px that may behave differently at finite p with 
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less or more accuracy. Though Eqs. (^) are equally satisfied for different choices of the 
coefficients Px the two states ipn and ip'^ may be much closer (less bias) for an optimal choice. 
The optimal choice that minimizes the distance {ipn ~ V'nl' ^^ fixed number p of correlation 
functions included in the SR, has not probably been found yet. We have attempted several 
choices for the reference weights Wj of Eq. ( P^D (with wj = Wj) but until now no 
significant improvement of the simplest FN ones [|l^] has been obtained. 



VIII. DETAILS OF THE ALGORITHM 

In this section the fiow chart of the GFMCSR algorithm is briefiy sketched. As described 
in the Appendix it is possible to work without the extra constant shift A and apply directly 
e"^"^, the usual imaginary time propagator, to filter out the ground state from the chosen 
trial wavefunction ipx. 

For practical purposes, the algorithm can be divided into three steps, 1) the Green 
function (OF) evolution, 2) the SR and 3) the measurements of physical mixed average 
correlation functions. These three steps are iterated until a satisfactory statistical accuracy 
is obtained for the latter quantities. 

The algorithm works with a finite number M of walkers. Starting from the first walker, 
corresponding conventionally to the index j = 1, the basic steps of the algorithm are de- 
scribed below: 

1. In the GF evolution, the exact propagator e"^^"^ and the FN one e'^" ^"^ are applied 
statistically for a given imaginary time interval At. In practice this can be done by 
setting initially Ati = At and repeating the following steps until Ati > 0: 

(a) Given the configuration of the walker, Xj, the quantities E^^, Vsf{xj) and H^^J^. 
Eqs. (P, p!6| , p!5| ) are evaluated. Then the interval At^ during which the walker is 
expected to perform only diagonal moves (see Appendix ^ is computed using 
the relation At^ = min(Arj, In^/TTf^), where ^ is a random number between and 
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1 and TTrf = limA->oo Alnp^ = E^j — H^^J^, according to Eq. (pi|). 

(b) Ati is updated Ar; — * Ar^— Ar^ and the walker weights {wj,Wj ) are muhiphed 
respectively by e(--^-.-(i+^)^-/(^^))^^'* and e'^^^^^^. Then if Ar^ > a new config- 
uration x'j 7^ Xj is chosen according to the probability table defined only by the 
normalized off-diagonal matrix elements oi px',x-, 

Px',Xj 



l^x'^Xj Px',Xj 

and the weight Wj is multiplied by s^'.^x (0)- The GF evolution then restarts 
from (a). Otherwise, if Ati = the GF evolution for the walker j terminates and 
the algorithm proceeds for the next walker starting from (1). 

2. After that all the walkers {wj,Wj ,Xj) have been propagated for the total imaginary 
time interval At the SR can be applied. The mixed averages O^. = {iPg\0\xj) / {ipclxj) 
are computed at the end of such propagation for the chosen set of operators O^. With 
these quantities both Ogjj = J2jWj 0^./J2Wj and the covariance matrix Sk,k' in 
Eq. ( P^ ) are evaluated. By using the latter quantities in the linear system (|C^), 
the coefficients a^ are easily computed and the table Px is determined according to 
Eq. ( |C^ . At this stage the reconfiguration procedure of the walkers can be eventually 
performed, i.e., the new M configurations of the walkers are chosen among the old 
ones according to the probability \Px\/J2k IPxJ- 

3. The mixed averages of the physical observables Oj and the quantity 

Hk'^k Hk \Pxk\ 

M EkPx, 
needed for the calculation of the statistical averages, are stored. The walker weights 
are set to Wj = sgnpx^ and Wj = 1, and the GF evolution can continue from step 
(1), starting again from the first walker. 

In the practical implementation of the algorithm the FN dynamic can be worked out at 
fixed 7, where 7 has to be a non-zero number otherwise the exact GF could not be sampled 
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(see Eqs. (|T^,|T5D). On the other hand for 7 = the FN is more accurate. A compromise 
is to work with 7 = 0.5 fixed. Another choice is to implement few runs with different 
non-zero 7 and try to extrapolate the results for 7 = 0, which should represent the most 
accurate calculation. Typically this extra effort is not necessary because there is a very weak 
dependence of the results upon 7. However the extrapolation to 7 — *> is an interesting 
possibility for the extension of the method to continuous models, since, in this case, there is 
no practical way to cross the nodes with a variational FN approach (as shown in Appendix 
for the lattice case). 

IX. THE LIMIT OF SMALL Ar AND LARGE NUMBER OF WALKERS 

In this section some general properties of the GFMCSR technique are discussed and 
explicitly tested on the J1 — J2 Heisenberg model 

^ = Ji ^ Si ■ S, + J2 E Si ■ S, , (44) 

(id) {{i,j)) 

where Sj are the s-1/2 operators sitting on the sites of a square lattice. Ji and J2 are 
the (positive) antiferromagnetic superexchange couplings between nearest and next-nearest- 
neighbors pairs of spins respectively. In the following we will consider finite square clusters 
of A^ sites with periodic boundary conditions. We use the same guiding wavefunction of 



Ref. |T0| and report here some test results useful to understand the crucial dependence of 
GFMCSR on the number of walkers M and the frequency of the SR Ar (the distance in 
imaginary time between two consecutive SR). In fact, after the selection of a given number 
p of correlation functions in Eqs. (|35D, the results depend only on the number of walkers M 
and the frequency of reconfiguration Ar. In the limit of large number of walkers, at fixed 
p, the algorithm has the important property that the fiuctuations of the coefficients ak and 
O'^ in Eq. (^) are obviously vanishing, because they depend on "averages" of a very large 
number of samples of many different walkers, implying that these fiuctuations are decreasing 
with 1/vM. In this limit it is possible to recover an important property of the FN: if the 
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guiding wavefunction is exact, the FN averages O^ are also exact. In fact suppose we begin 
to apply the propagator e~^'^ starting at r = from the exact samphng of the ground 
state ipo determined by FN with the exact guiding wavefunction ipc = V^o- Then at any 



Markov iteration n, before the SR is apphed, both the mixed average correlation functions 



calculated with the FN weights w^^^ {{O'') = ^ ^ ^fp ) and the weights with arbitrary 



. W ■ 



signs w {{O'') = -^ — —) sample statistically the true quantum average {iPq\0^\iIjo) . If, for 



large M, we can neglect statistical fluctuations of these averages, then by Eq. (|36| ) ak = 
and the SR algorithm just replace the weights Wj (with sign problem) with the FN weights 
Wj , which also sample tpo exactly if -ipc = ^o- This means that the SR approach does not 
affect this important property of the FN, at least in the limit M ^ oo. 

Another reason to work in the limit M — > cxd is the following. In this limit it is not 
necessary to put in the SR conditions (|36|) operators O^ that vanish for some symmetry that 
is satisfied both by the Hamiltonians H and the FN one H^^^ . In fact if the coefficients Px are 
defined in terms of operators O^ that conserve the mentioned symmetries (e.g. translation 
invariance, rotation by 90° degree of the lattice etc.) by definition Eqs. (^) are satisfied for 
all the remaining non-symmetric operators which have vanishing expectation value due to 
symmetry constraints (such as e.g. an operator that changes sign for a rotation operation 
which is a symmetry of the Hamiltonians H and H'^^^). In this case both sides of Eqs. (|35D 



are zero by such symmetry constraints. Moreover for M — ^ cx) the statistical fiuctuations are 
neglectable and for the same reason also Eqs. (|36|) are automatically satisfied with vanishing 
ak for the above mentioned non-symmetric operators. In this limit it is therefore useless to 
include non-symmetric operators in the SR (|36|) . 

Finally it is interesting that in this important limit M -^ oo, within the assumption 
that we can neglect the fiuctuations of a^ and O^^j, the SR depends only on the propagated 
states ipn'^^i^) a-^d ipnix). In fact given the state '?/'„(x) and the FN one ipn^^i^)^ then the 
state ipn{x) after the SR will be 
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(45) 

rn"'{x) = |<(X)| 

where now the ak are uniquely determined by the conditions (0), whereas the normahzation 
constant C = ^^ J}f , and, finally, ip^^-^' replace the FN propagated state ■^/'^-^-'^ after the 
SR (due to the condition Wj = \Wj\). In this hmit the dynamic described by the SR 
constraints is therefore perfectly defined and has a meaning, which can be computed even 
in an exact calculation without the Monte Carlo sampling. 

The way the computed results depend on the number of walkers is shown in Fig. |l|, as 
a function of the number of correcting factors. As it is evident for large number of walkers 
(M — ^ oo) the correcting factors do not play any role and the estimate with minimum 
statistical error is obtained by simply ignoring the correcting factors. This is actually a 
common approach in GFMC, to consider a large number of walkers so that the bias of 
the finite walker population becomes neglectable, and typically decreasing as 1/M (see e.g. 
Fig. 0). However from the picture it is also evident that for large enough M the predicted 
results obtained by including or by neglecting the correcting factors are both consistent. 
The convergence to the M — > oo limit is however faster for the first method. Thus the 
inclusion of the correcting factors G^ in Eq. (^O]), though increasing the error bars, may be 
useful to reach the M ^ oo limit with a smaller number of walkers. The fact that the two 
types of extrapolation to infinite M - the one including the correcting factors and the one 
neglecting them - converge to the same value (see Fig. Q) shows that the theoretical limit 
when (^5]) holds can be reached with a reasonable number of walkers, much smaller than 
the dimension of the Hilbert space. 

The other parameter that affects the accuracy of the SR approach is the imaginary time 
distance At between two consecutive SR. It is then natural to ask whether by increasing the 
frequency of the reconfigurations, one reaches a well defined dynamical limit for At -^ 0. 
This is important since, due to the sign problem for large size the time interval At has to be 
decreased at least by a factor inversely proportional to the system size, because the average 
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walker sign vanishes exponentially ~ e~^'''^ with an exponent A^ which diverges with the 
system size. Different calculations, performed for different sizes can be compared only when 
the finite At error (the difference between At —>■ and finite At) is neglectable. 

As shown in Fig. ^, whenever the simulation is stable for At -^ the limit At -^ can 
be reached with a linear extrapolation. This property can be easily understood since in the 
limit of large number of walkers the variation of the average correlation functions Eq.(^0|) 
both for the FN dynamic and the exact dynamic in a time interval between two consecutive 
SR differ clearly by 0(Ar). 

In order to show more clearly how the method is working and systematically correcting 
the FN we have implemented a slightly different but more straightforward "Release Nodes 
technique" ||15[. We first apply the standard FN (with 7 = 0, see Eq. (0)) for a given 



number of walkers M and for long simulation time. We store the M-walkers configurations, 
after some equilibration at time interval large enough to allow uncorrelated and independent 
samples of the FN ground state. In the second step we recover each of these M-walker 
configurations and apply GFMCSR for a fixed imaginary time r, so that we can see how the 
energy expectation value evolves from the FN to a more accurate determination. Typically 
one obtains a reasonable behavior for these curves that always coincides with the exact 
dynamic in the initial part where an exact sampling of the sign is possible. However for 
large imaginary time, and exceedingly small At and large number of walkers, some instability 
may occur leading to results clearly off, as shown in Fig. ^. In this case the reason of the 
instability is due to the fact that the correlation functions S'^(g) = -^Y,i,j SfS^e"^^"^'^^ 
which we have used in the SR (p = 9) [jlO[, introduce some uncontrolled fluctuations for the 
momentum Q = (vr, vr) relevant for the antiferromagnetic order parameter. If we include 
in the SR technique also the spin isotropic operator corresponding to the order parameter 
m''"^ = -^Y^ijSi ■ Sje'^'^^^~^^ and the total spin square (p = 11) this instability disappears 
(see Fig.^, stable results, not shown in the picture, are obtained even without the total spin 
square, i.e. with p = 10). This is a reasonable effect since the order parameter has important 
fluctuations in all spin directions. 
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X. CONCLUSIONS 

In this paper we have tried to describe in detail a recently proposed technique GFMCSR, 
that allows to work within a controlled accuracy with the ground state energy and with 
related mixed average correlation functions even for models where the conventional Quantum 
Monte Carlo technique cannot be used for the well known sign problem. 

This method is rather general, in principle convergence is achieved within an arbitrary 
accuracy if a sufficiently large number p of correlation functions is constraint to be equal be- 
fore and after the SR, the basic statistical step used to stabilize the sign problem instability. 
In order to minimize the number p of correlation functions used in the SR, one is limited 
to use an empirical approach, based on physical intuition, and/or by comparison with ex- 
act results obtained at finite size with the exact diagonalization technique. Typically the 
fundamental ingredient that we have found important for strongly correlated Hamiltonians 
is the "locality". The most useful correlation functions are the short range ones contained 
in the Hamiltonian H. A more successful example is the application of the method to the 
Heisenberg model on the triangular lattice [|16| where a remarkable accuracy is obtained 



by including also the short range correlation functions generated by the application of the 
square Hamiltonian. Here we report a table (see Table I) with all the values of the ground 
state energy per site, the total spin square and the antiferromagnetic order parameter m^"^ 
obtained with VMC, FN and GFMCSR (for two different p), up to A^ = 108. This method 
to increase systematically p, by including in the SR the short range correlation functions 
generated by H, H^ ■ ■ -, does not seem general enough. For instance it does not work for the 
Ji — J2 Heisenberg model where the inclusion of long range operators in the SR Eqs. (|35D 
such as the spin-spin correlation function S^Sj at large distance |i — j| is crucial to improve 
the accuracy of the method, whereas the short range ones generated by H^ do not give any 
significant improvement. 

Similarly to FN the GFMCSR is size- consistent (see Fig. |^). At fixed p a given accuracy 
is expected in the average correlation functions, accuracy which looks weakly dependent on 
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the system size and different from the variational guess even in the thermodynamic hmit. 
This is a very important property of the present algorithm because the stability of the 
average sign at fixed p allows a polynomial complexity of the algorithm as a function of the 
system size. The algorithm, however, is typically a large factor (~ 100) more expensive than 
the standard FN as far as the computational time is concerned, for a given statistical error 
on correlation functions. 

Until now the method has been extended rather successfully to several models: the men- 
tioned J1 — J2 and triangular lattice Heisenberg models, the t — J model Jl^ and preliminary 



results show that similar improvement of the standard FN can be obtained also for the 



Hubbard model [|1^]. In the latter case it is worth to mention that a different approach, the 
Constrained Path Monte Carlo |19[ (CPMC) represents also a very good remedy for the sign 
problem disease at least for intermediate coupling {U/t < 8). On the other hand different 
schemes to get rid of the "sign problem" for continuous systems were previously proposed 



and successfully applied to small electron systems. pO|] 

Although the GFMCSR is far from being the definite solution of the sign problem in 
the Monte Carlo simulation, it certainly represents an interesting possibility to alleviate 
this instability even for large system sizes. Its extension to continuous systems and also to 
CPMC is indeed straightforward, even though, in these cases, the possibility to cross the 
nodal surface in a variational way (see Appendix P) is not possible at present. 
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APPENDIX A: PROPERTY OF A STOCHASTIC MATRIX 

In this appendix we remind some properties of a stochastic matrix Px'^x- The stochastic 
matrices are square matrices that have all non- negative matrix elements px'^x and satisfy the 
normalization condition 

Ep-',- = i> (ai) 

x' 

for each column matrix index x. We assume also that the number of row and column 
indices are finite and that each index x is connected to any other x' by at least one sequence 
Px',xiPxi,x2 ' ' 'PxN,x of non-zero matrix elements of p. 

The stochastic matrices are generally non-symmetric and their eigenvalues may be also 
complex. For each eigenvalue there exist a left J2x' '4'l{x')Px',x = ^'^l{x) and a corresponding 
right eigenvector J2xPx',x'ipR{x) = \iPr{x'). A very simple left eigenvector is the constant 
one iPl{x) = 1, that by property (^T|) has eigenvalue A = 1. We will show in the following 
that this is actually the maximum eigenvalue because: i) to each right eigenvector iPr{x) of 
p corresponds an eigenvalue \, which is hounded by one |A| < 1 . 

In fact, be iPr{x) a generic (complex or real) right eigenvector of p 

X 

by taking the complex modulus of both sides of the previous equation and summing over x' 
we obtain 

x' x' X X x' X 

where in the above inequality we have interchanged the summation indices and used the 
elementary bound for the complex modulus | J2x^x\ < J2x \^x\ for arbitrary numbers Zx = 
Px',x4'r{x). This immediately gives: 

|A| < 1 . 

Obviously the equality sign holds if, for each x, \ J2x Zx\ = J2x kx|, which implies that given 
a right eigenvector with maximum eigenvalue A = 1, the real positive definite vector |'?/'ij(x)| 
is also a right eigenvector with maximum eigenvalue. 
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Now we will show that: ii) the maximum right eigenvector is unique. In fact suppose that 
there are two right eigenvectors ipi and 1^2 with A = 1, then by linearity also tpi—ailj2 is a right 
eigenvector with A = 1 and the complex constant a can be chosen to give tpi — a^2 = for a 
given index xq. On the other hand using the property derived previously also \i/ji{x)—aip2{x)\ 
is a right maximum eigenvector that vanishes for x = Xq. Using iteratively the definition of 
a right eigenvector 

Xl^^'.^lV^iW ~ «V^2(a;)| = |V^i(a;') - aip2ix')\ , 

X 

starting from x' = xq, we arrive easily to derive that for all the index x connected to xo by 
non-zero sequence of matrix elements Pxo,xiPxi,x2 ' ' 'Pxn,x 

|?/'i(x) -a'ilj2ix)\ = 0. 

Since by hypothesis all the possible indices are connected to xq by at least one such a 
sequence, we derive ipi = aip2, which means that ipi and 1IJ2 are the same eigenvector, which 
contradicts the initial hypothesis. Thus the maximum right eigenvector is unique. 

Collecting the above properties, the maximum right eigenvector iPr{x) of a stochastic 
matrix can be chosen real and positive. Then it is simple to show that the iteration of the 
stochastic matrix 

converges for large n to this maximum right eigenvector with an exponentially decreasing 
error oc 7", with 7 < 1 being the modulus of largest eigenvalue of p, different from the 
maximum one. 

APPENDIX B: PROOF OF THE UPPER BOUND 



Here we follow the paper |12] to prove rigorously the upper bound property of the ground 



state energy for H^^^ . We want to show that the prescription given in Eqs. (|n|,|T5|) for H^^^ 
leads to an upper bound for the ground state energy of H. When importance sampling is 
used it is important to change slightly the definition of the sign- flip term as in ([TBI) : 
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ipa{x')H^i .^/il)c{^)>^ and x'^x 

We now take any state 

X 

and we compare its energy with respect to H and to iJ^-^-^: 



(Bl) 



(B2) 



(B3) 



/\E can be written explicitly in terms of the matrix elements of H, using the definitions 



given in Eqs. (|14|,|15|,|BlD 



sf 



sf 



Vg[x) 



(B4) 



where the notation sf indicates conventionally the summation over the off-diagonal elements 
such that iI)g{x)Hx,x' / '4^g{x') > 0. In this double summation each pair of configurations x 
and x' occurs twice. We combine these terms and rewrite ( p4D as a summation over pairs: 



sf 



AE=(l+7) 5^ H^^^, 

{x,x') 



l^(x)|' ^|4^ + mx')f ^^ - ij{xyij{x') - ij{x'yi^{x) 



^g{x) 



i^oix') 



(B5) 



Denoting by sH{x,x') the sign of the matrix element H^.^./, and using the fact that for all 
terms in this summation the condition iPg{x')Hxi ^^ipcix) > is satisfied, we can finally write 

AE as 

2 



sf 



AE = (l + 7) 5] \Hx,x, 

{x,x') 



■^{x] 



\ 



tpcix') 



1pG{X) 



sH{x,x')ip{x'] 



\ 



iPg{x) 



iPg{x') 



(B6) 



Obviously, AE is positive for any wavefunction ip. Thus the ground state energy of H'^-I'f is 
an upper bound for the ground state energy of the original Hamiltonian H. 

Now the GFMC method can calculate the exact ground state energy Eq and wavefunc- 
tion ip'^-f^ of H^^^ , without any sign problem. Hence: -Eq > {il)^^^\H\il)^ff) > Eq, where the 
second inequality follows from the usual variational principle. We conclude therefore that 
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the FN energy is an upper bound to the true ground state energy. One can easily verify that 
{iPg\H\iPg) = {''Pg\H^^^\'4^g)^ and thus one can be sure that the GFMC procedure improves 
on the energy of the guiding wavef unction: E^^ < {^jJclH^-^^ipc) = {'4^g\H\iI)g)- 

Note that the standard "lattice FN" approach |T^ is obtained for the particular param- 
eter 7 = 0. 

APPENDIX C: PROOF OF EXISTENCE AND UNICITY OF SOLUTION FOR 

THE RECONFIGURATION 



In this appendix we prove that given the p + 1 SR conditions ( PE| ) the elements of the 
table Px are uniquely determined for each walker configuration (w, x). 
We define here the quantity 

v] = {Ol^ - O)) , (CI) 

for each configuration j, where O^ = ^ ^ p is the average value over the reference weights, 
Wj, of the operator considered, labeled by the number k. The reference weights Wj are 
restricted to be strictly positive but arbitrary functions of all the FN weights {Wj} the 
exact weights {wj} and the configurations {xj}. It is easy to show that, in order that 

P.,=w^^{l + J2a,v^) (C2) 

k 

allows to satisfy the SR conditions (p5D , it is sufficient that a^ are determined by the simple 
linear equation 

2^Sk,k'<yk' = -^ — - , (C3) 

k' Ej Wj 

where 

^— \ f k k' 

Sk,k' = J y (C4) 

is the covariance matrix of the operators O'' over the reference weights Wj . The solution to 
( P3|) is possible if the determinant of Sk^k' is non-vanishing. Since s represents an overlap 
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matrix defined with a non-singular scalar product {v'^\v^ ) = — L^^ ^j^ as w^ are positive, 
its determinant is always non-zero provided the vectors v'^ are linearly independent. Thus, 
in the latter case, the solution to ( |C3| ) exists and is unique. 

On the other hand suppose that among the p vectors v^ only p' <p are linearly indepen- 
dent. Thus the remaining p — p' vectors can be written as linear combination of p' linearly 
independent ones (henceforth we assume that these linearly independent vectors are labeled 
by the consecutive indices /c = 1, ■ ■ ■ ,p') 



V, 



p 

k' _ Sr^ ^k',k 



J 

fc=l 



E 4 ^ • , (C5) 



for k' > p', where x^ are suitable coefficients. The same previous considerations allow to 
satisfy the first p' SR conditions as for Eq. ( |C3| ) a unique solution exists if we restrict all the 
sums for k,k < p', and Px is determined only by the first p' linearly independent vectors in 
(P2D. With the determined p^ it is obvious that 



EjPx^_ Ej 



WiV 



^ (C6) 



is verified for k = 1,- ■ ■ ,p'. 

On the other hand we can easily show that all the remaining SR conditions (|C6D for 
k' > p' are identically satisfied. In fact, in this case the LHS of Eq. ([C^ ) can be manipulated 



as follows, using definition (|C5|) 



where in the intermediate steps we have used that 

for the first k < p' conditions. Thus the SR conditions determine uniquely Px in any case 
and this conclude the important statement of this Appendix. 
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1. Remark 

With the above definitions it is also possible to show that p^ remains unchanged for any 
linear transformation of the operator set. Namely, suppose we consider the new operators 

O'' = Y.Lk,,kO'' + Pk' (C8) 

k 

in the SR conditions, where the real matrix L is assumed to have no n- vanishing determinant. 
Within this assumption it is simple to show that p^ will remain unchanged. 

In fact the new set of operators will define a new covariance matrix between the new 
vectors 

vf = Y.L,,^,v^ , (C9) 

k 

i.e., V = Lv, s = LsL^ , where L^ is the transposed of L and the set of new equations 



1 )f^ 

j 



,^ ^ . _ Ej WjV 

}_. Sk,k'Oik' — —^ ' 

V Si ^3 

is obviously satisfied by 

a = {L-^fa , (CIO) 

where a is the solution of the SR conditions before the transformation (|C8|) . Whenever the 
number p' of linearly independent v^ is less than p, also the number of linearly independent 
v^ will be p' as L is non-singular. The solutions a and a, as described previously, refer 
therefore to the first p' components, and all the matrix involved, such as L and s are in this 
case restricted to this subspace. 

Then, by Eq. ( piO|) and Eq. ( |C9|) , it easily follows that the new coefficients px = 
Wj{l + Yjk'^kv'j) = wj(l + Yjk^^kV^) = Pxj, which finally proves the statement of this 
remark. 

APPENDIX D: THE LIMIT A ^ oo FOR THE POWER METHOD 

The constant A, which defines the the Green function Gx',x = ^Sx',x — H^'^x and the FN 
one C^^^ ([T3|) has to be taken large enough to determine that all the diagonal elements of 
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G^-^-^ are non- negative (by definition the off-diagonal ones of C^-^^ are always non- negative) . 
This requirement often determines a very large constant shift which increases with larger 
size and is not known a priori. The trouble in the simulation may be quite tedious, as if 
for the chosen A a negative diagonal element is found for C^-^^ , one needs to increase A and 
start again with a completely new simulation. The way out is to work with exceedingly 
large A, but this may slow down the efficiency of the algorithm as in the stochastic matrix 
Px\x the probability to remain in the same configuration pd may become very close to one 

Vd = ' ■ ^ ^ ^^ ^ ^ (Dl) 



where Vsf is given in Eq. ([TBD and E^ is the local energy Eq. (Jg) that do not depend on A 
given the configuration x. 

Following Ref. [0 the problem of working with large A can be easily solved with no loss 
of efficiency. We report this simple idea applied to our particular algorithm at fixed number 
of walkers. If A is large it is possible to take a large value of kp (of order A) iterations 
between two consecutive reconfigurations, because in most iterations the configuration x is 
not changed. The idea is that one can determine a priori, given pd what is the probability 
t{k) to make k diagonal moves before the first acceptance of a new configuration with x' ^ x. 
This is given by t{k) = Pd{^—Pd) for A; = 0, ■ ■ ■ , n^ — 1 and t(n/) = p^' if no off-diagonal moves 
are accepted during the ni trials that are left to complete the loop without reconfigurations. 

It is a simple exercise to show that, in order to sample t{k) one needs one random number 
< .^ < 1, so that the stochastic integer number k can be computed by the simple formula 

k = mmim, [p^]) , (D2) 

where the brackets indicate the integer part. During the kp iterations one can iteratively 
apply this formula by bookkeeping the number of iterations rii that are left to complete 
the loop without reconfigurations. At the first iteration ni = kp, then k is extracted using 
(P^), and the weights {w,w'^-^^) of the walker are updated according to k diagonal moves 
and ii k < rii a. new configuration is extracted randomly according to the off-diagonal matrix 
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elements oipx',x- The weights are correspondingly updated for this off-diagonal move, and 
finally, ii k < rii, rii is changed to ni — k — 1, so that one can continue to use Eq. (|D2| ) until 
all the kp steps are executed for each walker. 

The interesting thing of this method is that it can be readily generalized for A ^ oo by 
increasing kp with A, namely kp = [A At], where At represents now exactly the imaginary 
time difference between two consecutive reconfigurations where the exact propagator e~^^'^ 
or e"^" ^"^ is applied statistically. 
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TABLES 





N 


VMC 


FN 


SR(p = 2) 


SR{p = 7) 


Exact 


eo 


12 


-0.5981 


-0.6083(1) 


-0.6085(1) 


-0.6105(1) 


-0.6103 




36 


-0.5396 


-0.5469(1) 


-0.5534(1) 


-0.5581(1) 


-0.5604 




48 


-0.5366(1) 


-0.5426(1) 


-0.5495(1) 


-0.5541(1) 






108 


-0.5333(1) 


-0.5387(1) 


-0.5453(1) 


-0.5482(1) 




o2 
'-'tot 


12 


0.235 


0.0111(2) 


0.006(4) 


-0.002(4) 


0.00 




36 


1.71 


1.20(1) 


0.65(1) 


0.02(1) 


0.00 




48 


2.55(1) 


2.12(2) 


1.44(1) 


0.23(3) 


0.00 




108 


6.36(4) 


5.66(3) 


4.35(4) 


2.7(1) 


0.00 


mt2 


12 


0.9241 


0.9286(1) 


0.9210(2) 


0.9132(6) 


0.9109 




36 


0.7791 


0.7701(4) 


0.7659(2) 


0.7512(3) 


0.7394 




48 


0.7496(3) 


0.7243(5) 


0.7177(2) 


0.7080(5) 






108 


0.6338(7) 


0.6182(4) 


0.6040(3) 


0.5836(5) 





TABLE L Variational estimate (VMC) and mixed averages (FN, SR and Exact) of the ground 
energy per site, the total spin square and the order parameter for the triangular Heisenberg an- 
tiferromagnet for various system sizes. SR data are obtained using the short range correlation 



functions generated hy H {p = 2) and H^ {p = 7) reported in Ref. [16|. All the values reported 
in this table are obtained with large enough M and 1/Ar, practically converged in the limit of 
Ar ^ and M -^ infty. 
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FIGURES 
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40 60 

L 

FIG. 1. Dependence on the number L of correcting factors of the estimated ground state 
energy per site for A^ = 64 and J2 = 0.5 obtained with the GFMCSR technique (Ar = 0.01) with 
M = 200 (triangles), 1500 (squares) and 10000 (circles). 
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1/M 
FIG. 2. Ground state energy per site for J2 = 0.5 obtained for different clusters and different 

number of walkers. Empty dots are data obtained with zero correcting factors while full dots refer 

to converged values in L. 
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At 
FIG. 3. Dependence of the ground state energy per site on the imaginary time step Ar 

obtained for J2=0.5 and A^ = 36 with the GFMCSR technique by using in the SR the energy 

(full dots), all S^{q), the spin square and the order parameter m^'^ (empty dots). The number of 

walkers was fixed to M = 10000, so that the finite-M bias can be neglected on this scale. The 

lower horizontal axis coincides with the exact diagonalization result. 
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FIG. 4. Stable (upper curve) and unstable (lower curve) imaginary time evolution of the 
GFMCSR estimates of the ground energy per site for J2 = 0.5 and the A^ = 36 cluster. The 
horizontal line indicates the exact result. 
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0.001 0.002 0.003 0.004 0.005 

>f-3/2 

FIG. 5. Finite size scaling of the GS energy per site for J2 = 0.5 obtained with the FN and 



GFMCSR technique apphed reconfiguring the Hamiltonian (SRe). 
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